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Abstract 

The high computational cost involved in modeling of the progressive fracture simulations 
using large discrete lattice networks stems from the requirement to solve a new large set of 
linear equations every time a new lattice bond is broken. To address this problem, we pro- 
pose an algorithm that combines the multiple-rank sparse Cholesky downdating algorithm 
with the rank-p inverse updating algorithm based on the Sherman-Morrison- Woodbury for- 
mula for the simulation of progressive fracture in disordered quasi-brittle materials using 
discrete lattice networks. Using the present algorithm, the computational complexity of 
solving the new set of linear equations after breaking a bond reduces to the same order 
as that of a simple backsolve (forward elimination and backward substitution) using the 
already LU factored matrix. That is, the computational cost is 0(nnz(L)), where nnz(L) 
denotes the number of non-zeros of the Cholesky factorization L of the stiffness matrix A. 
This algorithm using the direct sparse solver is faster than the Fourier accelerated precondi- 
tioned conjugate gradient (PCG) iterative solvers, and eliminates the critical slowing down 
associated with the iterative solvers that is especially severe close to the critical points. 
Numerical results using random resistor networks substantiate the efficiency of the present 
algorithm. 
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1 Introduction 



Progressive damage evolution leading to failure of disordered quasi-brittle mate- 
rials has been studied extensively using various types of discrete lattice models 
. Large-scale numerical simulation of these lattice networks 
in which the damage is accumulated progressively by breaking one bond at a time 
until the lattice system falls apart has often been hampered due to the fact that a 
new large set of linear equations has to be solved everytime a lattice bond is bro- 
ken. Since the number of broken bonds at failure, n./, increases with increasing 
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lattice system sizes, L, i.e., rif ~ 0(L 17 ), numerical simulation of large lattice 
systems becomes prohibitively expensive. Furthermore, in fracture simulations us- 
ing discrete lattice networks, ensemble averaging of numerical results is necessary 
to obtain a realistic representation of the lattice system response. This further in- 
creases the computational cost associated with modeling fracture simulations in 
disordered quasi-brittle materials using large discrete lattice networks. 

Fourier accelerated PCG iterative solvers Jjl [HI II 1 have been used in the past 
for simulating the material breakdown using large lattices. However, these methods 
do not completely eliminate the critical slowing down associated with the iterative 
solvers close to the critical point. As the lattice system gets closer to macroscopic 
fracture, the condition number of the system of linear equations increases, thereby 
increasing the number of iterations required to attain a fixed accuracy. This becomes 
particularly significant for large lattices. Furthermore, the Fourier acceleration tech- 
nique is not effective when fracture simulation is performed using central-force and 
bond-bending lattice models [Holl . 



This study presents an algorithm that combines the multiple-rank sparse Cholesky 
downdating scheme with the rank-p inverse updating scheme of the stiffness matrix, 
which effectively reduces the computational bottleneck involved in re-solving the 
new set of equations after everytime a bond is broken. In this paper, we consider a 
random threshold model problem, where a lattice consists of fuses having the same 
conductance, but the bond breaking thresholds, i c , are based on a broad (uniform) 
probability distribution, which is constant between and 1. This relatively simple 
model has been extensively used in the literature fH I3> IE SI H> El 0] for simulat- 
ing the fracture and progressive damage evolution in brittle materials, and provides 
a meaningful benchmark for comparing different algorithms. A broad thresholds 
distribution represents large disorder and exhibits diffusive damage leading to pro- 
gressive localization, whereas a very narrow thresholds distribution exhibits brittle 
failure in which a single crack propagation causes material failure. Periodic bound- 
ary conditions are imposed in the horizontal direction to simulate an infinite system 
and a constant voltage difference (displacement) is applied between the top and the 
bottom of lattice system. The simulation is initiated with a triangular lattice of in- 
tact fuses of size L x L, in which disorder is introduced through random breaking 
thresholds. The voltage V across the lattice system is increased until a fuse (bond 
breaking) burns out. The burning of a fuse occurs whenever the electrical current 
(stress) in the fuse (bond) exceeds the breaking threshold current (stress) value 
of the fuse. The current is redistributed instantaneously after a fuse is burnt. The 
voltage is then gradually increased until a second fuse is burnt, and the process is 
repeated. Each time a fuse is removed, the electrical current is redistributed and 
hence it is necessary to re-solve Kirchhoff equations to determine the current flow- 
ing in the remaining bonds of the lattice. This step is essential for determining the 
fuse that is going to burn up under the redistributed currents. Therefore, numeri- 
cal simulations leading to final breaking of lattice system network are very time 
consuming especially with increasing lattice system size. 
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1.1 Summary of the Proposed Algorithm 

The algorithm presented in this paper reduces the computational complexity of 
obtaining the solution x n , after the n th bond is broken, to a backsolve using the 
already existent factorization of the stiffness matrix A m , and p = (n — m) vector 
updates. The algorithm is based on the well known Shermon-Morrison-Woodbury 
formula 1 11211 for obtaining the inverse of the new stiffness matrix A"^ (after the 
(n + l) th fuse is burnt) from the old stiffness matrix inverse A" 1 through a rank-one 
update. Infact, the algorithm is such that if the inverse of the lattice stiffness at any 
stage (m = 0, 1, 2 ■ ■ ■) of analysis A" 1 is available, then all subsequent analysis 
involving (n = m + l,m + 2, ■ ■ ■) burnt fuses can be carried out using p — (n — m) 
vector updates. However, since most often the inverse of the stiffness matrix is 
rarely ever explicitly calculated, the algorithm additionally requires a backsolve 
using the already existent factored matrix A m . The backsolve operation is further 
simplified by the fact that it is performed on a trivial load vector and hence the 
solution can be obtained easily. 

Based on the above description of the algorithm presented in this paper, given the 
factorization of the matrix A m , the computational cost involved in all the subse- 
quent steps (n = m + 1, m + 2, • • •) is a backsolve using the already factored 
matrix, and p — (n — m) vector updates. The computational complexity of the 
backsolve is 0(nnz(L m )), where nnz(L m ) denotes the number of non-zeros of 
the Cholesky factorization L m of A m . The computational complexity of p vector 
updates is 0(p n do f), where n do f denotes the number of degrees of freedom in the 
system. As p increases, it is possible that the computational cost associated with 
the p vector updates exceeds the cost involved in the factorization of the matrix 
A n . Under these circumstances, it is advantageous to obtain the factorization L n of 
the new stiffness matrix A n , and use this L n for all the subsequent backsolve analy- 
sis steps, until the computational cost associated with the vector updates once again 
exceeds the stiffness factorization cost. Using the algorithm presented in the paper, 
it is not necessary to re-factorize the new stiffness matrix A n . Instead, we adopt 
the multiple-rank update of the sparse Cholesky factorization algorithm llii \\M 
for updating the L m — > L„. This multiple-rank update of L m to obtain the new 
factorization L n is computationally cheaper compared with the direct factorization 
L n of the new stiffness matrix A n HHEI. 



2 Proposed Algorithm 

In the following, we describe the updating scheme for the inverse of the stiffness 
matrix in the case of scalar random fuse model after a fuse has been burnt. A similar 
procedure can be applied for central-force and beam models [15]. 
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Let A n represent the stiffness matrix of the random fuse network system in which 
n number of fuses are either missing (random dilution) or have been burnt during 
the analysis. Let us also assume that a fuse ij (the (n + l) th fuse) is burnt when the 
externally applied voltage is increased gradually. In the above description, % and j 
refer to the global degrees of freedom connected by the fuse before it is broken. For 
the scalar random fuse model, the degrees of freedom i and j are also equivalent to 
the node i and node j connected by the fuse before it is broken. The new stiffness 
matrix A n+1 of the lattice system after the fuse ij is burnt is given by 

A n+ i = A n - kij vv* (1) 



where 



v* = | • • • 1 1 • • • } (2) 

and is the conductance of the fuse ij before it is broken. After breaking the fuse 
ij, the electrical current in the network is redistributed instantaneously. The redis- 
tributed current values in the network are calculated by re-solving the Kirchhoff 
equations, i.e., by solving the new set of equations formed by the matrix A n+1 . 
This procedure is very time consuming since a new set of equations (inverse of 
A n+ i for n — 0, 1, 2, • • •) need to be solved everytime after breaking the (n + l) th 
fuse. However, significant computational advantages can be gained if the inverse 
of A n+ i is obtained simply by updating the inverse of A n . This is achieved by us- 
ing the well known Shermon-Morrison- Woodbury formula for inverting the rank-p 
update of a matrix. Thus, the inverse A~^ x of Eq. (1) can be expressed as 



A^ 1 - 



A„ + ki 



VlVl 



(1 - k i:j v*u)_ 



(3) 



where 



u = A n 1 v = A n l {i _ 3) 



■th 



X 



j th ) columns of A^ 1 ] (4) 



Hence, the inverse of the stiffness matrix of the lattice system after breaking the 
(n + l) th fuse ij is obtained simply by a rank-one update of the inverse of the 
stiffness matrix before the fuse is broken. Further, if the inverse of the matrix A n 
is available explicitly, then the vector u can be obtained trivially from the i th and 
j th columns of A" 1 . In particular, this implies that if the inverse of the matrix A n 
is available explicitly at any stage n — 0, 1, 2 • • • of analysis, then the redistributed 
currents in all subsequent stages of analysis involving m — n + 1, n + 2, • • • burnt 
fuses can be obtained in a trivial fashion from the column vectors of A" 1 and the 
vectors u p , where p — 1, 2, • • • , (m — n). However, since the inverse of the stiffness 
matrix A n is not usually calculated explicitly, the vector u is obtained using the 
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already factorized A n matrix through a backsolve operation (forward reduction 
and backward substitution) on the vector v (Eq. (4)). 

REMARK 1: Without loss of generality, when the fuse that is broken is attached to 
a constrained/prescribed degree of freedom j, the vector v is given by 



{o-i-o} 



and 



(5) 



i th columns of A: 



(6) 



In the case of periodic boundary conditions, consider the case of a broken fuse jk 
that is attached to a slave degree of freedom k whose master degree of freedom is 
i. Under these circumstances, the methodology presented earlier is applicable in a 
straightforward manner if it is understood that breaking the fuse jk is equivalent to 
breaking the fuse ij. 

REMARK 2: The load vector b n+i will differ from the load vector b n only if the 
(n + l) th broken fuse ij is attached to a prescribed degree of freedom, where a 
constant voltage difference is imposed. Once again, for presentation purposes, let 
us assume that j is such a prescribed degree of freedom. Then the load vector b n+1 
is given by 

b„+i = K + w (7) 
where 

w* = k i;j j 1 • • • j (8) 

If neither i nor j is a prescribed degree of freedom, then w = 0. 

Before breaking the (n + l) th fuse, the solution vector x n is obtained by solving 
the Kirchhoff equations 

A n x n = b n (9) 



After breaking the (n + l) th fuse that connects the i th and j th degrees of freedom, 
the updated solution vector x n+i is obtained by solving the new set of Kirchhoff 
equations 

A n+ ix n+1 = b n+1 (10) 
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Substituting Eqs. (3,7) and (9) into the solution of Eq. (10) and simplifying the 
result, we have 



X„+l — A^jbn+i 



A n 1 + kij 



(b n + w) 



(1 - v*u)_ 

x„ + 0u (11) 



where 



[3 = a (u*b n+ i) — kij if i or j is prescribed 
— a ^u*b„+i^ otherwise (12) 

and 



<x= jz k :° i , (B) 

(1 - hj v*u) 



The only unknown in Eq. (11) is the column vector u, which can be obtained 
through a backsolve operation using either Eq. (4) or Eq. (6). Furthermore, it is 
not necessary to explicitly assemble the matrix A n and perform factorization to do 
the backsolve operation. Instead, we can use the already factorized matrix A rn to 
obtain the vector u. In the above description, m < n and denotes the latest broken 
bond at which the factorization L m of A m is available. To see this clearly, let us 
first decompose the matrix A" 1 into A" 1 and a matrix C such that 

A; 1 = A- 1 + C (14) 
where 



p=(n-m) t 

C = E *t n I t \ (15) 
ti (1 - h vjui) 

Due to the amount of the storage requirement (~ 0(n 2 do ^)), and the computational 
cost associated in evaluating the Eq. (15) (~ 0(n 2 do j)), the matrix C is never ex- 
plicitly calculated or stored. Instead, the vectors u; for I = 1, 2, • • • , (n — m) are 
stored, and the (j th — i th ) column of C is evaluated as 



p=(n-m) / \ 

C «->= g ^ (l4viu,) "' <16) 
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where and refer to the i th and j th components of the vector u^. Equation (16) 
reduces the storage and computational cost to (~ 0(p ridof)) and (~ 0(p rido})) 
operations, respectively. Even with this modification, the storage and computational 
requirements can become prohibitively expensive as the number of updates, p, in- 
creases, and hence it is necessary to limit the maximum number of vector updates 
between two successive factorizations to a certain maxupd. That is, it is necessary 
to perform or update the factorization of the stiffness matrix A at regular intervals. 

Instead of re-factorizing the stiffness matrix A after every maxupd steps, it is more 
effective to update the factorization L m using the multiple-rank sparse Cholesky 



factorization update algorithm I13L 11411 . This multiple-rank update of L m to obtain 
the new factorization L n+1 , after breaking the (n + l) th fuse, is computationally 
cheaper compared with the direct factorization L n+1 of the new stiffness matrix 
A n+ i J13IH41I . We use the multiple-rank downdate algorithm presented in [13, 14] 
to obtain the new Cholesky factorization L n +] f rom the existing Cholesky factor 
L m . The multiple-rank downdate algorithm H13L is based on the analysis and 
manipulation of the underlying graph structure of the stiffness matrix A and on 
the methodology presented in Gill et al. 1 16, 17J] for modifying a dense Cholesky 
factorization. This algorithm incorporates the change in the sparsity pattern of L 
and is optimal in the sense that the computational time required is proportional 
to the number of changing non-zero entries in L. In particular, since the breaking 
of fuses is equivalent to removing the edges in the underlying graph structure of 
stiffness matrix A, the new sparsity pattern of the modified L must be a subset of 
the sparsity pattern of the original L. Denoting the sparsity pattern of L by C, we 
have 



C m D C n V m < n (17) 
Therefore, we can even use the modified dense Cholesky factorization update (al- 



gorithm 5 in the reference Davis et al. tUP ) and work only on the non-zero entries 
in L. Furthermore, since the changing non-zero entries in L depend on the i th and 
jth degrees f freedom of the fuse ij that is broken, it is only necessary to modify 
the non-zero elements of a submatrix of L. 

The multiple-rank update of the sparse Cholesky factorization is computationally 
superior to an equivalent series of rank-one updates since the multiple-rank update 
makes one pass through L in computing the new entries, while a series of rank-one 
updates require multiple passes through L [ 14]. The multiple-rank update algorithm 
updates the Cholesky factorization L m of the matrix A m to L n+ i of the new matrix 
A n+ i, where A rt+1 = A m + aYY*, and Y represents a n do f x p rank-p matrix. 
The computational cost involved in breaking the (n + l) th fuse ij is simply a back- 
solve operation (0(nnz(L m ))) on a load vector given by Eq. (2) using the already 
factored matrix A m , [n + 2 — m) vector updates, and one vector inner product. 

The optimum number of steps between successive factorizations of the matrix A is 
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determined by minimizing the computatioal cpu time required for the entire analy- 
sis. Let tf ac and t upd denote the average cpu time required for performing/updating 
the factorization A m and the average cpu time required for a single rank-1 update of 
the solution U( n+1 _ m ), respectively. Note that the evaluation of U( n+ i_ m ) requires 
(n — m) vector updates. Let the estimated number of steps for the lattice system 
failure be n sieps . Then, the total cpu time required for solving the linear system of 
equations until the lattice system failure is given by 



^upd 

itfac 

1 {jl steps Tlfac) steps 



nfactfac + ~ — t upd (18) 



where nj ac denotes the number of factorization until lattice system failure. The 
optimum number of factorizations, n opi __/ ac , for the entire analysis is obtained by 
minimizing the function ty. The maximum number of vector updates, maxupd, 
between successive factorizations is estimated as 

7 (Tlsteps Tlopt— fac) m 

maxupd = — - — (19) 

^opt—fac 



3 Numerical Simulation Results 



In the following, we consider two alternate forms of the algorithm presented in this 
paper. These two solver types are 

• Solver Type A: Given the factorization L m of A m , we use rank-1 sparse Cholesky 



update/downdate [13] to update the factorization L n+1 (0(nnz(L n )) for all sub- 
sequent values ofn = m,m + l,'". Once the factorization L n+1 of A n+1 is ob- 
tained, the solution vector x n+ i is obtained by a backsolve operation (0(nnz(L n+ i)). 
• Solver Type B: Given the factorization L m of A m , the algorithm evaluates the 
new solution vector x n+1 , after the [n+l) th fuse is burnt, using Eq. (11) {0{nnz{L m ) 
+ (n + 2 — m) vector updates). Instead of refactorizing the matrix after maxupd 
steps, we use rank-p sparse Cholesky update/downdate 111411 to obtain the factor- 
ization L m+maxupd of the matrix A m+maxupd (0(nnz(L m )). 

The above two alg orithms are benchmarked against the PCG iterative solvers, in 
which optimal 12111 circulant matrices are used as preconditioners to the 



Laplacian operator (Kirchhoff equations). The Fourier accelerated PCG presented 



in 



1 (All lH is not optimal in the sense described in [ 18, 19, 20, 21], and hence it is 
expected to take more number of CG iterations compared with the optimal circulant 
preconditioners. 
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In the numerical simulations using solver types A and B, the supernodal Cholesky 
factorization is performed using the TAUCS solver library (http://www.tau.ac. il^stoledo/taucs). 
In these simulations, the maximum number of vector updates, maxupd, is cho- 
sen to be a constant for a given lattice size. We choose maxupd = 25 for L = 
{4, 8, 16, 24, 32}, maxupd = 50 for L = 64, and maxupd = 100 for L = 
{128,256,512}. For L = 512, maxupd is limited to 100 due to memory con- 
straints. By keeping the maxupd value constant, it is possible to realistically com- 
pare the computational cost associated with different solver types. Moreover, the 
relative cpu times taken by these algorithms remains the same even when the sim- 
ulations are performed on different platforms. 

Tables 1 and 2 present the cpu and wall-clock times taken for one configuration 
(simulation) using the solver types A and B, respectively. These tables also indi- 
cate the number of configurations, N con f ig , over which ensemble averaging of the 
numerical results is performed. The cpu and wall-clock times taken by the opti- 
mal circulant matrix preconditioned iterative solver is presented in Tables 3. For 
iterative solvers, the number of iterations presented in Tables 3 denote the average 
number of total iterations taken to break one intact lattice configuration until it falls 
apart. 

Based on the results presented in Tables 1-3, it is clear that for modeling the break- 
down of disordered media as in starting with an intact lattice and successive break- 
ing of bonds until the lattice system falls apart, the solver types A and B based 
on direct solvers are superior to the Fourier accelerated iterative PCG solver tech- 
niques. It should be noted that for larger lattice systems, limitations on the available 
memory of the processor may decrease the allowable maxupd value, as in the case 
of L = 512 using solver type B. However, this is not a concern for simulations 
performed using solver type A. 

Using the solver type A, we have performed numerical simulations on two-dimensional 
triangular and diamond (square lattice inclined at 45 degrees between the bus bars) 
lattice networks. Table 4 presents the number of broken bonds at peak load, n p , 
and at fracture, rif, for each of the lattice sizes considered. In addition, Table 4 also 
presents the number of configurations, N con f ig , over which statistical averaging is 
performed for different lattice sizes. The numerical results presented in Tables 1-3 
are performed on a single processor of Cheetah (27 Regatta nodes with thirty two 
1.3 GHz Power4 processors each), the eighth fastest supercomputer in the world 
(http://www.ccs.ornl.gov). However, the numerical simulation results presented in 
Table 4 are performed on Eagle (184 nodes with four 375 MHz Power3-II pro- 
cessors) supercomputer at the Oak Ridge National Laboratory to run simulations 
simultaneously on more number of processors. Figure 1 presents the snapshots of 
progressive damage evolution for the case of a broadly distributed random thresh- 
olds model problem in a triangular lattice system of size L = 512. 
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4 Conclusions 



The paper presents an algorithm based on rank-one update of the inverse of the stiff- 
ness matrix and the multiple-rank downdating of the sparse Cholesky factorization 
for simulating fracture and damage evolution in disordered quasi-brittle materials 
using discrete lattice networks. Using the proposed algorithm, the average compu- 
tational cost associated with breaking a bond reduces to the same order as that of a 
simple backsolve (forward elimination and backward substitution) operation using 
the already LU factored matrix. This algorithm based on direct solver techniques 
eliminates critical slowing down observed in fracture simulations using the con- 
ventional iterative schemes. Numerical simulations on random resistor networks 
demonstrate that the present algorithm is computationally superior to the com- 
monly used Fourier accelerated preconditioned conjugate gradient iterative solver. 

For analysis of fracture simulations using discrete lattice networks, ensemble av- 
eraging of numerical results is necessary to obtain a realistic representation of the 
lattice system response. In this regard, for very large lattice systems with large 
number of system of equations, this methodology is especially advantageous as 
the LU factorization of the system of equations can be performed using a parallel 
implementation on multiple processors. Subsequently, this factored LU decompo- 
sition can then be distributed to each of the processors to continue with independent 
fracture simulations that only require less intensive backsolve operations. 

Acknowledgment 

This research is sponsored by the Mathematical, Information and Computational 
Sciences Division, Office of Advanced Scientific Computing Research, U.S. De- 
partment of Energy under contract number DE-AC05-00OR22725 with UT-Battelle, 
LLC. The first author wishes to thank Ed F. D'Azevedo for many helpful discus- 
sions and excellent suggestions. 



10 



References 

[I] L. de Arcangelis, S. Redner, and H. J. Herrmann. A random fuse model for 
breaking processes. Journal of Physics (Paris) Letters, 46(13):585-590, 1985. 

[2] M. Sahimi and J. D. Goddard. Elastic percolation models for cohesive me- 
chanical failure in heterogeneous systems. Physical Review B, 33:7848-7851, 
1986. 

[3] P. M. Duxbury, P. D. Beale, and P. L. Leath. Size effects of electrical break- 
down in quenched random media. Physical Review Letters, 57(8): 1052-1055, 
1986. 

[4] P. M. Duxbury, P. L. Leath, and P. D. Beale. Breakdown properties of 
quenched random systems: The random-fuse network. Physical Review B, 
36:367-380, 1987. 

[5] A. Hansen and S. Roux. Statistical toolbox for damage and fracture, pages 
17-101. Springer, New York, 2000. 

[6] H. J. Herrmann and S. Roux (eds). Statistical Models for the Fracture of 
Disordered Media. North-Holland, Amsterdam, 1990. 

[7] M. Sahimi. Non-linear and non-local transport processes in heterogeneous 
media from long-range correlation percolation to fracture and materials break- 
down. Physics Reports, 306:213-395, 1998. 

[8] Bikas K. Chakrabarti and L. Gilles Benguigui. Statistical Physics of Fracture 
and Breakdown in Disordered Systems. Oxford Science Publications, Oxford, 
1997. 

[9] G. G. Batrouni, A. Hansen, and M. Nelkin. Fourier acceleration of relaxation 
processes in disordered systems. Physical Review Letters, 57:1336-1339, 
1986. 

[10] G. G. Batrouni and A. Hansen. Fourier acceleration of iterative processes in 
disordered-systems. Journal of Statistical Physics, 52:747-773, 1988. 

[I I] G. G. Batrouni and A. Hansen. Fracture in three-dimensional fuse networks. 
Physical Review Letters, 80:325-328, 1998. 

[12] G. H. Golub and C. F. van Loan. Matrix Computations. The Johns Hopkins 
University Press, 1996. 

[13] Timothy A. Davis and William W. Hager. Modifying a sparse Cholesky fac- 
torization. SIAMJ. Matrix Anal. AppL, 20(3):606-627, 1999. 

[14] Timothy A. Davis and William W. Hager. Multiple-rank modifications of a 
sparse Cholesky factorization. SI AM J. Matrix Anal. AppL, 22(4):997-1013, 
2001. 

[15] Phani Kumar V. V. Nukala and Srdan Simunovic. unpublished. 

[16] P. E. Gill, G. H. Golub, W. Murray, and M. A. Saunders. Methods for modi- 
fying matrix factorizations. Math. Comp., 28:505-535, 1974. 

[17] P. E. Gill, W. Murray, and M. A. Saunders. Methods for computing and mod- 
ifying the LDV factors of a matrix. Math. Comp., 29:1051-1077, 1975. 

[18] T. Chan. An optimal circulant preconditioner for Toeplitz systems. SIAM J. 
Sci. Stat. Comput., 9:766-771, 1988. 

[19] Raymond H. Chan. Circulant preconditioners for Hermitian Toeplitz systems. 



11 



SIAMJ. Matrix Anal. AppL, 10:542-550, 1989. 

[20] R. Chan and T. Chan. Circulant preconditioners for Elliptic problems. Nu- 
merical Linear Algebra Applications, 1:77-101, 1992. 

[21] Raymond H. Chan and Michael K. Ng. Conjugate gradient methods for 
Toeplitz systems. SIAM Review, 38(3):427-482, 1996. 



12 



Table 1 

Computational cost associated with solver type A 



Size 


CPU(sec) 


Wall(sec) 


Simulations 


32 


0.592 


0.687 


20000 


64 


10.72 


11.26 


4000 


128 


212.2 


214.9 


800 


256 


5647 


5662 


96 


512 


93779 


96515 


16 



Table 2 

Computational cost associated with solver type B 



Size 


CPU(sec) 


Wall(sec) 


Simulations 


32 


0.543 


0.633 


20000 


64 


11.15 


12.01 


4000 


128 


211.5 


214.1 


800 


256 


6413 


6701 


96 



Table 3 

Computational cost associated with optimal circulant PCG 



Size 


CPU(sec) 


Wall(sec) 


Iterations 


Simulations 


32 


11.66 


12.26 


25469 


20000 


64 


173.6 


178.8 


120570 


1600 


128 


7473 


7725 


622140 


128 
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Table 4 

Number of broken bonds at peak and at failure 



L 


config 


time 
(seconds) 


Triangular 


Diamond 


¥ 


J 


¥ 


J 


4 


50000 


0.002 


13 


19 


9 


14 


8 


50000 


0.006 


41 


54 


26 


37 


16 


50000 


0.042 


134 


168 


80 


107 


24 


50000 


0.186 


276 


335 


161 


208 


32 


50000 


0.592 


465 


554 


268 


337 


64 


50000 


10.72 


1662 


1911 


942 


1126 


128 


12000 


212.2 


6068 


6766 


3406 


3901 


256 


1200 


5647 


22572 


24474 


12571 


13846 


512 


200 


93779 


84487 


89595 
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Fig. 1. Snapshots of damage in a typical triangular lattice system of size L = 512. Number 
of broken bonds at the peak load and at failure are 83995 and 89100, respectively, (a)-(i) 
represent the snapshots of damage after breaking n b number of bonds, (a) n b = 25000 (b) 
n b = 50000 (c) n b = 75000 (d) n b = 80000 (e) n b = 83995 (peak load) (f) n b = 86000 
(g) n b = 87000 (h) n b = 88000 (i) n b = 89100 (failure) 
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